Create a new analysis directory...
[1] FALSE
[1] FALSE
[1] FALSE
[1] FALSE
[1] FALSE
[1] FALSE
[1] "/Users/slaan3/git/CirculatoryHealth/AE_20211201_YAW_SWVANDERLAAN_HDAC9"
[1] "_archived" "1. AE_20211201_YAW_SWVANDERLAAN_HDAC9.nb.html"
[3] "1. AE_20211201_YAW_SWVANDERLAAN_HDAC9.Rmd" "1. AEDB.CEA.baseline.nb.html"
[5] "1. AEDB.CEA.baseline.Rmd" "2. SNP_analyses.nb.html"
[7] "2. SNP_analyses.Rmd" "20220319.HDAC9.AEDB.CEA.baseline.RData"
[9] "20220319.HDAC9.AESCRNA.results.RData" "20220319.HDAC9.bulkRNAseq.additional_figures.RData"
[11] "20220319.HDAC9.bulkRNAseq.main_analysis.RData" "20220319.HDAC9.bulkRNAseq.preparation.RData"
[13] "20230301.HDAC9.bulkRNAseq.additional_figures.RData" "20230511.HDAC9.bulkRNAseq.additional_figures.RData"
[15] "20230531.HDAC9.bulkRNAseq.additional_figures.RData" "3.1 bulkRNAseq.preparation.nb.html"
[17] "3.1 bulkRNAseq.preparation.Rmd" "3.2 bulkRNAseq.main_analysis.nb.html"
[19] "3.2 bulkRNAseq.main_analysis.Rmd" "3.3 bulkRNAseq.additional_figures.nb.html"
[21] "3.3 bulkRNAseq.additional_figures.Rmd" "4. scRNAseq.nb.html"
[23] "4. scRNAseq.Rmd" "AE_20211201_YAW_SWVANDERLAAN_HDAC9.Rproj"
[25] "AnalysisPlan" "HDAC9"
[27] "images" "LICENSE"
[29] "README.html" "README.md"
[31] "references.bib" "renv"
[33] "renv.lock" "scripts"
[35] "SNP"
source(paste0(PROJECT_loc, "/scripts/functions.R"))install.packages.auto("readr")
install.packages.auto("optparse")
install.packages.auto("tools")
install.packages.auto("dplyr")
install.packages.auto("tidyr")
install.packages.auto("naniar")
# To get 'data.table' with 'fwrite' to be able to directly write gzipped-files
# Ref: https://stackoverflow.com/questions/42788401/is-possible-to-use-fwrite-from-data-table-with-gzfile
# install.packages("data.table", repos = "https://Rdatatable.gitlab.io/data.table")
library(data.table)
install.packages.auto("tidyverse")
install.packages.auto("knitr")
install.packages.auto("DT")
install.packages.auto("MASS")
# install.packages.auto("Seurat") # latest version
# Install the devtools package from Hadley Wickham
install.packages.auto('devtools')
install.packages.auto("haven")
install.packages.auto("sjlabelled")
install.packages.auto("sjPlot")
install.packages.auto("labelled")
install.packages.auto("tableone")
install.packages.auto("ggpubr")This notebook contains additional figures of the project “Plaque expression levels of HDAC9 in association with plaque vulnerability traits and secondary vascular events in patients undergoing carotid endarterectomy: an analysis in the Athero-EXPRESS Biobank.”.
# load(paste0(PROJECT_loc, "/",Today,".",PROJECTNAME,".bulkRNAseq.main_analysis.RData"))
load(paste0(PROJECT_loc, "/20220319.",PROJECTNAME,".bulkRNAseq.main_analysis.RData"))We need to get the ‘conventional unit’ versions of cholesterols.
AERNASE.clin.hdac9 <- merge(AERNASE.clin.hdac9,
subset(AEDB.CEA, select = c("STUDY_NUMBER",
"risk614",
"LDL_finalCU", "HDL_finalCU", "TC_finalCU", "TG_finalCU")),
by.x = "STUDY_NUMBER", by.y = "STUDY_NUMBER", sort = TRUE, all.x = TRUE)We want to create per-age-group figures median ± interquartile range.
# ?ggpubr::ggboxplot()
compare_means(HDAC9 ~ Gender, data = AERNASE.clin.hdac9, method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9,
x = c("Gender"),
y = "HDAC9",
xlab = "gender",
ylab = "HDAC9 (normalized expression)",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Gender.pdf"), plot = last_plot())Saving 12 x 8 in image
library(dplyr)
AERNASE.clin.hdac9 <- AERNASE.clin.hdac9 %>% dplyr::mutate(AgeGroup = factor(case_when(Age < 55 ~ "<55",
Age >= 55 & Age <= 64 ~ "55-64",
Age >= 65 & Age <= 74 ~ "65-74",
Age >= 75 & Age <= 84 ~ "75-84",
Age >= 85 ~ "85+")))
AERNASE.clin.hdac9 <- AERNASE.clin.hdac9 %>% dplyr::mutate(AgeGroupSex = factor(case_when(Age < 55 & Gender == "male" ~ "<55 males" ,
Age >= 55 & Age <= 64 & Gender == "male"~ "55-64 males",
Age >= 65 & Age <= 74 & Gender == "male"~ "65-74 males",
Age >= 75 & Age <= 84 & Gender == "male"~ "75-84 males",
Age >= 85 & Gender == "male"~ "85+ males",
Age < 55 & Gender == "female" ~ "<55 females" ,
Age >= 55 & Age <= 64 & Gender == "female"~ "55-64 females ",
Age >= 65 & Age <= 74 & Gender == "female"~ "65-74 females",
Age >= 75 & Age <= 84 & Gender == "female"~ "75-84 females",
Age >= 85 & Gender == "female"~ "85+ females")))
table(AERNASE.clin.hdac9$AgeGroup, AERNASE.clin.hdac9$Gender)
female male
<55 11 27
55-64 43 124
65-74 58 191
75-84 37 119
85+ 4 9
table(AERNASE.clin.hdac9$AgeGroupSex)
<55 females <55 males 55-64 females 55-64 males 65-74 females 65-74 males 75-84 females 75-84 males 85+ females
11 27 43 124 58 191 37 119 4
85+ males
9
Now we can draw some graphs of plaque target(s) levels per sex and age group as median ± interquartile range.
# ?ggpubr::ggboxplot()
compare_means(HDAC9 ~ AgeGroup, data = AERNASE.clin.hdac9, method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9,
x = c("AgeGroup"),
y = "HDAC9",
xlab = "Age groups (years)",
ylab = "HDAC9 (normalized expression)",
color = "AgeGroup",
palette = "npg",
# add = "median_iqr")
add = c("median_iqr", "jitter")) +
stat_compare_means(aes(group = AgeGroup), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.AgeGroup.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ AgeGroup, group.by = "Gender", data = AERNASE.clin.hdac9, method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9,
x = c("AgeGroup"),
y = "HDAC9",
xlab = "Age groups (years) per gender",
ylab = "HDAC9 (normalized expression",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
# add = "median_iqr")
add = c("median_iqr", "jitter")) +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.AgeGroup_perGender.pdf"), plot = last_plot())Saving 12 x 8 in image
We want to create figures of target(s) levels stratified by hypertension/blood pressure, and use of anti-hypertensive drugs.
library(dplyr)
AERNASE.clin.hdac9 <- AERNASE.clin.hdac9 %>% mutate(SBPGroup = factor(case_when(systolic < 120 ~ "<120",
systolic >= 120 & systolic <= 139 ~ "120-139",
systolic >= 140 & systolic <= 159 ~ "140-159",
systolic >= 160 ~ "160+")))
table(AERNASE.clin.hdac9$SBPGroup, AERNASE.clin.hdac9$Gender)
female male
<120 7 22
120-139 30 81
140-159 36 120
160+ 62 169
Now we can draw some graphs of plaque target(s) levels per sex and hypertension/blood pressure group as median ± interquartile range.
compare_means(HDAC9 ~ SBPGroup, data = AERNASE.clin.hdac9 %>% filter(!is.na(SBPGroup)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(SBPGroup)),
x = c("SBPGroup"),
y = "HDAC9",
xlab = "Systolic blood pressure (mmHg)",
ylab = "HDAC9 (normalized expression)",
color = "SBPGroup",
palette = "npg",
add = "jitter") +
stat_compare_means(aes(group = SBPGroup), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.SBPGroup.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ Hypertension.selfreport, data = AERNASE.clin.hdac9 %>% filter(!is.na(Hypertension.selfreport)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Hypertension.selfreport)),
x = c("Hypertension.selfreport"),
y = "HDAC9",
xlab = "Self-reported hypertension",
ylab = "HDAC9 (normalized expression)",
color = "Hypertension.selfreport",
palette = "npg",
add = "jitter") +
stat_compare_means(aes(group = Hypertension.selfreport), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Hypertension.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ Hypertension.drugs, data = AERNASE.clin.hdac9 %>% filter(!is.na(Hypertension.drugs)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Hypertension.drugs)),
x = c("Hypertension.drugs"),
y = "HDAC9",
xlab = "Hypertension medication use",
ylab = "HDAC9 (normalized expression)",
color = "Hypertension.drugs",
palette = "npg",
add = "jitter") +
stat_compare_means(aes(group = Hypertension.drugs), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.HypertensionDrugs.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ SBPGroup, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(SBPGroup)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(SBPGroup)),
x = c("SBPGroup"),
y = "HDAC9",
xlab = "Systolic blood pressure (mmHg) per gender",
ylab = "HDAC9 (normalized expression)",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.SBPGroup_byGender.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ Hypertension.selfreport, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(Hypertension.selfreport)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Hypertension.selfreport)),
x = c("Hypertension.selfreport"),
y = "HDAC9",
xlab = "Self-reported hypertension per gender",
ylab = "HDAC9 (normalized expression)",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Hypertension_byGender.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ Hypertension.drugs, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(Hypertension.drugs)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Hypertension.drugs)),
x = c("Hypertension.drugs"),
y = "HDAC9",
xlab = "Hypertension medication use per gender",
ylab = "HDAC9 (normalized expression)",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Hypertension.drugs_byGender.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ SBPGroup, group.by = "Hypertension.drugs", data = AERNASE.clin.hdac9 %>% filter(!is.na(SBPGroup) & !is.na(Hypertension.drugs)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(SBPGroup) & !is.na(Hypertension.drugs)),
x = c("SBPGroup"),
y = "HDAC9",
xlab = "Systolic blood pressure (mmHg) by medication use",
ylab = "HDAC9 (normalized expression)",
color = "Hypertension.drugs",
palette = c("#49A01D", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Hypertension.drugs), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.SBPGroup_byHypertensionDrugs.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ Hypertension.selfreport, group.by = "Hypertension.drugs", data = AERNASE.clin.hdac9 %>% filter(!is.na(Hypertension.selfreport) & !is.na(Hypertension.drugs)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Hypertension.selfreport) & !is.na(Hypertension.drugs)),
x = c("Hypertension.selfreport"),
y = "HDAC9",
xlab = "Self-reported hypertension per medication use",
ylab = "HDAC9 (normalized expression)",
color = "Hypertension.drugs",
palette = c("#49A01D", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Hypertension.drugs), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Hypertension.selfreport_byHypertensionDrugs.pdf"), plot = last_plot())Saving 12 x 8 in image
We want to create figures of target(s) levels stratified by hypercholesterolemia/LDL-levels, and use of lipid-lowering drugs.
risk614) group (no, yes)library(dplyr)
AERNASE.clin.hdac9 <- AERNASE.clin.hdac9 %>% mutate(LDLGroup = factor(case_when(LDL_finalCU < 100 ~ "<100",
LDL_finalCU >= 100 & LDL_finalCU <= 129 ~ "100-129",
LDL_finalCU >= 130 & LDL_finalCU <= 159 ~ "130-159",
LDL_finalCU >= 160 & LDL_finalCU <= 189 ~ "160-189",
LDL_finalCU >= 190 ~ "190+")))
table(AERNASE.clin.hdac9$LDLGroup, AERNASE.clin.hdac9$Gender)
female male
<100 45 142
100-129 25 73
130-159 18 42
160-189 9 21
190+ 2 9
require(sjlabelled)
AERNASE.clin.hdac9$risk614 <- to_factor(AERNASE.clin.hdac9$risk614)
# Fix plaquephenotypes
attach(AERNASE.clin.hdac9)
AERNASE.clin.hdac9[,"Hypercholesterolemia"] <- NA
AERNASE.clin.hdac9$Hypercholesterolemia[risk614 == "missing value"] <- NA
AERNASE.clin.hdac9$Hypercholesterolemia[risk614 == -999] <- NA
AERNASE.clin.hdac9$Hypercholesterolemia[risk614 == 0] <- "no"
AERNASE.clin.hdac9$Hypercholesterolemia[risk614 == 1] <- "yes"
detach(AERNASE.clin.hdac9)
table(AERNASE.clin.hdac9$risk614, AERNASE.clin.hdac9$Hypercholesterolemia)< table of extent 3 x 0 >
# AEDB.temp <- subset(AEDB, select = c("STUDY_NUMBER", "UPID", "Age", "Gender", "Hospital", "Artery_summary", "risk614", "Hypercholesterolemia"))
# require(labelled)
# AEDB.temp$Gender <- to_factor(AEDB.temp$Gender)
# AEDB.temp$Hospital <- to_factor(AEDB.temp$Hospital)
# AEDB.temp$Artery_summary <- to_factor(AEDB.temp$Artery_summary)
#
# DT::datatable(AEDB.temp[1:10,], caption = "Excerpt of the whole AEDB.", rownames = FALSE)
#
# rm(AEDB.temp)Now we can draw some graphs of plaque target(s) levels per sex and hypercholesterolemia/LDL-levels group, as well as stratified by lipid-lowering drugs users as median ± interquartile range.
compare_means(HDAC9 ~ LDLGroup, data = AERNASE.clin.hdac9 %>% filter(!is.na(LDLGroup)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(LDLGroup)),
x = c("LDLGroup"),
y = "HDAC9",
xlab = "LDL (mg/dL) per gender",
ylab = "HDAC9 (normalized expression))",
color = "LDLGroup",
palette = "npg",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.LDLGroups.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ LDLGroup, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(LDLGroup)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(LDLGroup)),
x = c("LDLGroup"),
y = "HDAC9",
xlab = "LDL (mg/dL) per gender",
ylab = "HDAC9 (normalized expression))",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.LDLGroups_byGender.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ Hypercholesterolemia, data = AERNASE.clin.hdac9 %>% filter(!is.na(Hypercholesterolemia)), method = "kruskal.test")Error in kruskal.test.default(x = mf[[1L]], g = mf[[2L]]) :
all observations are in the same group
ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Hypercholesterolemia)),
x = c("Hypercholesterolemia"),
y = "HDAC9",
xlab = "Diagnosed hypercholesterolemia",
ylab = "HDAC9 (normalized expression))",
color = "Hypercholesterolemia",
palette = "npg",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Hypercholesterolemia.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ Hypercholesterolemia, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(Hypercholesterolemia)), method = "kruskal.test")Error in names(grouped.d$data) <- .names :
'names' attribute [1] must be the same length as the vector [0]
ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Hypercholesterolemia)),
x = c("Hypercholesterolemia"),
y = "HDAC9",
xlab = "Diagnosed hypercholesterolemia per gender",
ylab = "HDAC9 (normalized expression))",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Hypercholesterolemia_byGender.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ Med.Statin.LLD, data = AERNASE.clin.hdac9 %>% filter(!is.na(Med.Statin.LLD)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Med.Statin.LLD)),
x = c("Med.Statin.LLD"),
y = "HDAC9",
xlab = "Lipid-lowering drug use",
ylab = "HDAC9 (normalized expression))",
color = "Med.Statin.LLD",
palette = "npg",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Med.Statin.LLD.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ Med.Statin.LLD, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(Med.Statin.LLD)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Med.Statin.LLD)),
x = c("Med.Statin.LLD"),
y = "HDAC9",
xlab = "Lipid-lowering drug use per gender",
ylab = "HDAC9 (normalized expression))",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Med.Statin.LLD_byGender.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ LDLGroup, group.by = "Med.Statin.LLD", data = AERNASE.clin.hdac9 %>% filter(!is.na(LDLGroup) & !is.na(Med.Statin.LLD)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(LDLGroup) & !is.na(Med.Statin.LLD)),
x = c("LDLGroup"),
y = "HDAC9",
xlab = "LDL (mg/dL) per LLD use",
ylab = "HDAC9 (normalized expression))",
color = "Med.Statin.LLD",
palette = c("#49A01D", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Med.Statin.LLD), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.LDLGroups_byMed.Statin.LLD.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ Hypercholesterolemia, group.by = "Med.Statin.LLD", data = AERNASE.clin.hdac9 %>% filter(!is.na(Hypercholesterolemia) & !is.na(Med.Statin.LLD)), method = "kruskal.test")Error in names(grouped.d$data) <- .names :
'names' attribute [1] must be the same length as the vector [0]
ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Hypercholesterolemia) & !is.na(Med.Statin.LLD)),
x = c("Hypercholesterolemia"),
y = "HDAC9",
xlab = "Diagnosed hypercholesterolemia per LLD use",
ylab = "HDAC9 (normalized expression))",
color = "Med.Statin.LLD",
palette = c("#49A01D", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Med.Statin.LLD), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.LDLGroups_byMed.Statin.LLD.pdf"), plot = last_plot())Saving 12 x 8 in image
We want to create figures of target(s) levels stratified by kidney function.
library(dplyr)
AERNASE.clin.hdac9 <- AERNASE.clin.hdac9 %>% mutate(eGFRGroup = factor(case_when(GFR_MDRD < 15 ~ "<15",
GFR_MDRD >= 15 & GFR_MDRD <= 29 ~ "15-29",
GFR_MDRD >= 30 & GFR_MDRD <= 59 ~ "30-59",
GFR_MDRD >= 60 & GFR_MDRD <= 89 ~ "60-89",
GFR_MDRD >= 90 ~ "90+")))
table(AERNASE.clin.hdac9$eGFRGroup, AERNASE.clin.hdac9$Gender)
female male
15-29 2 6
30-59 38 101
60-89 84 249
90+ 25 86
table(AERNASE.clin.hdac9$eGFRGroup, AERNASE.clin.hdac9$KDOQI)
No data available/missing Normal kidney function CKD 2 (Mild) CKD 3 (Moderate) CKD 4 (Severe) CKD 5 (Failure)
15-29 0 0 0 0 8 0
30-59 0 0 0 139 0 0
60-89 0 0 333 0 0 0
90+ 0 111 0 0 0 0
Now we can draw some graphs of plaque target(s) levels per sex and kidney function group as median ± interquartile range.
# Global test
compare_means(HDAC9 ~ eGFRGroup, data = AERNASE.clin.hdac9 %>% filter(!is.na(eGFRGroup)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(eGFRGroup)),
x = c("eGFRGroup"),
y = "HDAC9",
xlab = "eGFR (mL/min per 1.73 m2)",
ylab = "HDAC9 (normalized expression)",
color = "eGFRGroup",
palette = "npg",
add = "jitter") +
stat_compare_means(method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.EGFR.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ eGFRGroup, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(eGFRGroup)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(eGFRGroup)),
x = c("eGFRGroup"),
y = "HDAC9",
xlab = "eGFR (mL/min per 1.73 m2) per gender",
ylab = "HDAC9 (normalized expression)",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.EGFR_byGender.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ KDOQI, data = AERNASE.clin.hdac9 %>% filter(!is.na(KDOQI)), method = "kruskal.test")p1 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(KDOQI)),
x = c("KDOQI"),
y = "HDAC9",
xlab = "Kidney function (KDOQI)",
ylab = "HDAC9 (normalized expression)",
color = "KDOQI",
palette = "npg",
add = "jitter") +
stat_compare_means(aes(group = KDOQI), label = "p.format", method = "kruskal.test")
ggpar(p1 + rotate_x_text(45), legend = "right") rm(p1)
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.KDOQI.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ KDOQI, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(KDOQI)), method = "kruskal.test")p1 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(KDOQI)),
x = c("KDOQI"),
y = "HDAC9",
xlab = "Kidney function (KDOQI) per gender",
ylab = "HDAC9 (normalized expression)",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")
ggpar(p1 + rotate_x_text(45), legend = "right") rm(p1)
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.KDOQI_byGender.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ eGFRGroup, data = AERNASE.clin.hdac9 %>% filter(!is.na(eGFRGroup) & !is.na(KDOQI)), method = "kruskal.test")p1 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(eGFRGroup) & !is.na(KDOQI)),
x = c("eGFRGroup"),
y = "HDAC9",
xlab = "eGFR (mL/min per 1.73 m2) by KDOQI group",
ylab = "HDAC9 (normalized expression)",
color = "KDOQI",
palette = "npg",
add = "jitter") +
stat_compare_means(method = "kruskal.test")
ggpar(p1, legend = "right")rm(p1)
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.EGFR_KDOQI.pdf"), plot = last_plot())Saving 12 x 8 in image
We want to create figures of target(s) levels stratified by BMI.
library(dplyr)
AERNASE.clin.hdac9 <- AERNASE.clin.hdac9 %>% mutate(BMIGroup = factor(case_when(BMI < 18.5 ~ "<18.5",
BMI >= 18.5 & BMI < 25 ~ "18.5-24",
BMI >= 25 & BMI < 30 ~ "25-29",
BMI >= 30 & BMI < 35 ~ "30-35",
BMI >= 35 ~ "35+")))
# require(labelled)
# AERNASE.clin.hdac9$BMI_US <- as_factor(AERNASE.clin.hdac9$BMI_US)
# AERNASE.clin.hdac9$BMI_WHO <- as_factor(AERNASE.clin.hdac9$BMI_WHO)
# table(AERNASE.clin.hdac9$BMI_WHO, AERNASE.clin.hdac9$BMI_US)
table(AERNASE.clin.hdac9$BMIGroup, AERNASE.clin.hdac9$Gender)
female male
<18.5 3 2
18.5-24 46 162
25-29 68 220
30-35 18 55
35+ 6 12
table(AERNASE.clin.hdac9$BMIGroup, AERNASE.clin.hdac9$BMI_WHO)
No data available/missing Underweight Normal Overweight Obese
<18.5 0 5 0 0 0
18.5-24 0 0 208 0 0
25-29 0 0 0 287 0
30-35 0 0 0 0 73
35+ 0 0 0 0 18
Now we can draw some graphs of plaque MCP1 levels per sex and age group as median ± interquartile range.
# Global test
compare_means(HDAC9 ~ BMIGroup, data = AERNASE.clin.hdac9 %>% filter(!is.na(BMIGroup)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(BMIGroup)),
x = c("BMIGroup"),
y = "HDAC9",
xlab = "BMI groups (kg/m2)",
ylab = "HDAC9 (normalized expression)",
# color = "Gender",
# palette = c("#D5267B", "#1290D9"),
color = "BMIGroup",
palette = "npg",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.BMI.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ BMIGroup, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(BMIGroup)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(BMIGroup)),
x = c("BMIGroup"),
y = "HDAC9",
xlab = "BMI groups (kg/m2) per gender",
ylab = "HDAC9 (normalized expression)",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.BMI_byGender.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ BMIGroup, data = AERNASE.clin.hdac9 %>% filter(!is.na(BMIGroup) & !is.na(BMI_WHO)), method = "kruskal.test")p1 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(BMIGroup) & !is.na(BMI_WHO)),
x = c("BMIGroup"),
y = "HDAC9",
xlab = "BMI groups (kg/m2) per WHO categories",
ylab = "HDAC9 (normalized expression)",
color = "BMI_WHO",
palette = "npg",
add = "jitter") +
stat_compare_means(method = "kruskal.test")
ggpar(p1, legend = "right")rm(p1)
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.BMI_byWHO.pdf"), plot = last_plot())Saving 12 x 8 in image
We want to create figures of target(s) levels stratified by type 2 diabetes.
Now we can draw some graphs of plaque target(s) levels per sex and age group as median ± interquartile range.
# Global test
compare_means(HDAC9 ~ DiabetesStatus, data = AERNASE.clin.hdac9 %>% filter(!is.na(DiabetesStatus)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(DiabetesStatus)),
x = c("DiabetesStatus"),
y = "HDAC9",
xlab = "Diabetes status",
ylab = "HDAC9 (normalized expression)",
# color = "Gender",
# palette = c("#D5267B", "#1290D9"),
color = "DiabetesStatus",
palette = "npg",
add = c("median_iqr", "jitter")) +
stat_compare_means(label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Diabetes.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ DiabetesStatus, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(DiabetesStatus)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(DiabetesStatus)),
x = c("DiabetesStatus"),
y = "HDAC9",
xlab = "Diabetes status per gender",
ylab = "HDAC9 (normalized expression)",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = c("median_iqr", "jitter")) +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Diabetes_byGender.pdf"), plot = last_plot())Saving 12 x 8 in image
We want to create figures of target(s) levels stratified by smoking.
Now we can draw some graphs of plaque target(s) levels per sex and age group as median ± interquartile range.
# Global test
compare_means(HDAC9 ~ SmokerStatus, data = AERNASE.clin.hdac9 %>% filter(!is.na(SmokerStatus)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(SmokerStatus)),
x = c("SmokerStatus"),
y = "HDAC9",
xlab = "Smoker status",
ylab = "HDAC9 (normalized expression)",
# color = "Gender",
# palette = c("#D5267B", "#1290D9"),
color = "SmokerStatus",
palette = "npg",
add = c("median_iqr", "jitter")) +
stat_compare_means(label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Smoking.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ SmokerStatus, group.by ="Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(SmokerStatus)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(SmokerStatus)),
x = c("SmokerStatus"),
y = "HDAC9",
xlab = "Smoker status per gender",
ylab = "HDAC9 (normalized expression)",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = c("median_iqr", "jitter")) +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Smoking_byGender.pdf"), plot = last_plot())Saving 12 x 8 in image
We want to create figures of target(s) levels stratified by stenosis grade.
library(dplyr)
AERNASE.clin.hdac9 <- AERNASE.clin.hdac9 %>% mutate(StenoticGroup = factor(case_when(stenose == "0-49%" ~ "<70",
stenose == "0-49%" ~ "<70",
stenose == "50-70%" ~ "<70",
stenose == "70-90%" ~ "70-89",
stenose == "50-99%" ~ "90+",
stenose == "70-99%" ~ "90+",
stenose == "100% (Occlusion)" ~ "90+",
stenose == "90-99%" ~ "90+")))
table(AERNASE.clin.hdac9$StenoticGroup, AERNASE.clin.hdac9$Gender)
female male
<70 6 34
70-89 72 199
90+ 69 221
table(AERNASE.clin.hdac9$stenose, AERNASE.clin.hdac9$StenoticGroup)
<70 70-89 90+
missing 0 0 0
0-49% 2 0 0
50-70% 38 0 0
70-90% 0 271 0
90-99% 0 0 284
100% (Occlusion) 0 0 5
NA 0 0 0
50-99% 0 0 1
70-99% 0 0 0
99 0 0 0
Now we can draw some graphs of plaque target(s) levels per sex and age group as median ± interquartile range.
# Global test
compare_means(HDAC9 ~ StenoticGroup, data = AERNASE.clin.hdac9 %>% filter(!is.na(StenoticGroup)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(StenoticGroup)),
x = c("StenoticGroup"),
y = "HDAC9",
xlab = "Stenotic grade",
ylab = "HDAC9 (normalized expression)",
# color = "Gender",
# palette = c("#D5267B", "#1290D9"),
color = "StenoticGroup",
palette = "npg",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Stenosis.pdf"), plot = last_plot())Saving 12 x 8 in image
compare_means(HDAC9 ~ StenoticGroup, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(StenoticGroup)), method = "kruskal.test")ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(StenoticGroup)),
x = c("StenoticGroup"),
y = "HDAC9",
xlab = "Stenotic grade per gender",
ylab = "HDAC9 (normalized expression)",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.Stenosis_byGender.pdf"), plot = last_plot())Saving 12 x 8 in image
We want to create per-symptom figures.
library(dplyr)
table(AERNASE.clin.hdac9$AgeGroup, AERNASE.clin.hdac9$AsymptSympt2G)
Asymptomatic Symptomatic
<55 10 28
55-64 24 143
65-74 30 219
75-84 15 141
85+ 1 12
table(AERNASE.clin.hdac9$Gender, AERNASE.clin.hdac9$AsymptSympt2G)
Asymptomatic Symptomatic
female 15 138
male 65 405
table(AERNASE.clin.hdac9$AsymptSympt2G)
Asymptomatic Symptomatic
80 543
Now we can draw some graphs of plaque target(s) levels per symptom group as median ± interquartile range.
# ?ggpubr::ggboxplot()
my_comparisons <- list(c("Asymptomatic", "Symptomatic"))
compare_means(HDAC9 ~ AsymptSympt2G, data = AERNASE.clin.hdac9, method = "kruskal.test")p1 <- ggpubr::ggboxplot(AERNASE.clin.hdac9,
x = "AsymptSympt2G", y = "HDAC9",
title = "HDAC9 (normalized expression) levels per symptom",
xlab = "Symptoms",
ylab = "HDAC9 (normalized expression)",
color = "AsymptSympt2G",
# palette = c(uithof_color[16], uithof_color[23]),
palette = "npg",
add = "dotplot", # Add dotplot
add.params = list(binwidth = 0.1, dotsize = 0.3)
) +
stat_compare_means(comparisons = my_comparisons, method = "wilcox.test")
ggpar(p1, legend = c("right"), legend.title = "Symptoms")
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.AsymptSympt2G.pdf"), plot = last_plot())Saving 12 x 8 in image
rm(p1)compare_means(HDAC9 ~ AsymptSympt2G, group.by = "Gender", data = AERNASE.clin.hdac9, method = "kruskal.test")p1 <- ggpubr::ggboxplot(AERNASE.clin.hdac9,
x = "AsymptSympt2G", y = "HDAC9",
title = "HDAC9 (normalized expression) levels per symptom by gender",
xlab = "Symptoms",
ylab = "HDAC9 (normalized expression)",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "dotplot", # Add dotplot
add.params = list(binwidth = 0.1, dotsize = 0.3)
) +
stat_compare_means(aes(group = Gender), label = "p.format", method = "wilcox.test")
ggpar(p1, legend = c("right"), legend.title = "Symptoms")
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.AsymptSympt2G.byGender.pdf"), plot = last_plot())Saving 12 x 8 in image
rm(p1)# ?ggpubr::ggboxplot()
# compare_means(LRCH1 ~ eGFRGroup, data = AERNASE.clin %>% filter(!is.na(eGFRGroup)), method = "kruskal.test")
# ggpubr::ggboxplot(AERNASE.clin %>% filter(!is.na(eGFRGroup))
cat("Get summary statistics for target:\n")Get summary statistics for target:
summary(AERNASE.clin.hdac9$HDAC9) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 11.00 21.00 28.32 37.00 449.00
cat("\nCount number of values > 100 for target:\n")
Count number of values > 100 for target:
sum(AERNASE.clin.hdac9$HDAC9 > 100)[1] 12
cat("\nSetting values > 100 for target to 100.\n")
Setting values > 100 for target to 100.
temp <- AERNASE.clin.hdac9 %>%
filter(!is.na(AsymptSympt2G))
temp$HDAC9[temp$HDAC9 > 100] <- 100
cat("Get summary statistics for target after fixing values:\n")Get summary statistics for target after fixing values:
summary(temp$HDAC9) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 11.00 21.00 26.71 37.00 100.00
my_comparisons <- list(c("Asymptomatic", "Symptomatic"))
compare_means(HDAC9 ~ AsymptSympt2G, data = temp, method = "kruskal.test")
p1 <- ggpubr::ggboxplot(temp,
x = "AsymptSympt2G", y = "HDAC9",
title = "HDAC9 (normalized expression) levels per symptom",
xlab = "Symptoms",
ylab = "HDAC9 (normalized expression)",
color = "AsymptSympt2G",
# palette = c(uithof_color[16], uithof_color[23]),
palette = "npg",
add = "jitter", # Add dotplot
add.params = list(binwidth = 0.1, dotsize = 0.3)
) +
stat_compare_means(comparisons = my_comparisons, method = "wilcox.test")
ggpar(p1, legend = c("right"), legend.title = "Symptoms")
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.",TRAIT_OF_INTEREST,".HDAC9.plaque.AsymptSympt2G.noNA_limitCount100.pdf"), plot = last_plot())Saving 12 x 8 in image
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.",TRAIT_OF_INTEREST,".HDAC9.plaque.AsymptSympt2G.noNA_limitCount100.ps"), plot = last_plot())Saving 12 x 8 in image
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.",TRAIT_OF_INTEREST,".HDAC9.plaque.AsymptSympt2G.noNA_limitCount100.png"), plot = last_plot())Saving 12 x 8 in image
rm(p1, temp)We would also like to visualize the multivariable analyses results.
library(ggplot2)
library(openxlsx)
model1_target <- read.xlsx(paste0(OUT_loc, "/", Today, ".AERNASE.clin.hdac9.Bin.Uni.",TRAIT_OF_INTEREST,".RANK.Symptoms.MODEL1.xlsx"))
model2_target <- read.xlsx(paste0(OUT_loc, "/", Today, ".AERNASE.clin.hdac9.Bin.Multi.",TRAIT_OF_INTEREST,".RANK.Symptoms.MODEL2.xlsx"))
model1_target$model <- "univariate"
model2_target$model <- "multivariate"
models_target <- rbind(model1_target, model2_target)
models_targetForest plots.
dat <- data.frame(group = factor(c("Age, sex-adjusted", "Age, sex, and adjusted for risk factors"),
levels=c("Age, sex, and adjusted for risk factors", "Age, sex-adjusted")),
cen = c(models_target$OR[models_target$Predictor=="HDAC9"]),
low = c(models_target$low95CI[models_target$Predictor=="HDAC9"]),
high = c(models_target$up95CI[models_target$Predictor=="HDAC9"]))
fp <- ggplot(data = dat, aes(x = group, y = cen, ymin = low, ymax = high)) +
geom_pointrange(linetype = 2, size = 1, colour = c("#1290D9", "#49A01D")) +
geom_hline(yintercept = 1, lty = 2) + # add a dotted line at x=1 after flip
coord_flip(ylim = c(0.8, 1.7)) + # flip coordinates (puts labels on y axis)
xlab("Model") + ylab("OR (95% CI) for symptomatic plaques") +
ggtitle("Normalized HDAC9 expression") +
theme_minimal() # use a white background
print(fp)
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,".plaque.forest.pdf"), plot = fp)
rm(fp)We will plot the correlations of other cytokine plaque levels to the MCP1 plaque levels. These include:
In addition we will look at three metalloproteinases which were measured using an activity assay.
The proteins were measured using FACS and LUMINEX. Given the different platforms used (FACS vs. LUMINEX), we will inverse rank-normalize these variables as well to scale them to the same scale as the target(s)` plaque levels.
We will set the measurements that yielded ‘0’ to NA, as it is unlikely that any protein ever has exactly 0 copies. The ‘0’ yielded during the experiment are due to the limits of the detection.
# fix names
names(AEDB.CEA)[names(AEDB.CEA) == "VEFGA"] <- "VEGFA"
# fix names
names(AERNASE.clin.hdac9)[names(AERNASE.clin.hdac9) == "IL6"] <- "IL6rna"
cytokines <- c("IL2", "IL4", "IL5", "IL6", "IL8", "IL9", "IL10", "IL12", "IL13", "IL21",
"INFG", "TNFA", "MIF", "MCP1", "MIP1a", "RANTES", "MIG", "IP10", "Eotaxin1",
"TARC", "PARC", "MDC", "OPG", "sICAM1", "VEGFA", "TGFB")
metalloproteinases <- c("MMP2", "MMP8", "MMP9")
AERNASE.clin.hdac9 <- merge(AERNASE.clin.hdac9,
subset(AEDB.CEA, select = c("STUDY_NUMBER",
cytokines,
metalloproteinases)),
by.x = "STUDY_NUMBER", by.y = "STUDY_NUMBER", sort = TRUE, all.x = TRUE)
proteins_of_interest <- c(cytokines, metalloproteinases)
proteins_of_interest_rank = unlist(lapply(proteins_of_interest, paste0, "_rank"))
# make variables numerics()
AERNASE.clin.hdac9 <- AERNASE.clin.hdac9 %>%
mutate_each(funs(as.numeric), proteins_of_interest)
for(PROTEIN in 1:length(proteins_of_interest)){
# UCORBIOGSAqc$Z <- NULL
var.temp.rank = proteins_of_interest_rank[PROTEIN]
var.temp = proteins_of_interest[PROTEIN]
cat(paste0("\nSelecting ", var.temp, " and standardising: ", var.temp.rank,".\n"))
cat(paste0("* changing ", var.temp, " to numeric.\n"))
# AERNASE.clin.hdac9 <- AERNASE.clin.hdac9 %>% mutate(AERNASE.clin.hdac9[,var.temp] == replace(AERNASE.clin.hdac9[,var.temp], AERNASE.clin.hdac9[,var.temp]==0, NA))
AERNASE.clin.hdac9[,var.temp][AERNASE.clin.hdac9[,var.temp]==0.000000]=NA
cat(paste0("* standardising ", var.temp,
" (mean: ",round(mean(!is.na(AERNASE.clin.hdac9[,var.temp])), digits = 6),
", n = ",sum(!is.na(AERNASE.clin.hdac9[,var.temp])),").\n"))
AERNASE.clin.hdac9 <- AERNASE.clin.hdac9 %>%
mutate_at(vars(var.temp),
# list(Z = ~ (AERNASE.clin.hdac9[,var.temp] - mean(AERNASE.clin.hdac9[,var.temp], na.rm = TRUE))/sd(AERNASE.clin.hdac9[,var.temp], na.rm = TRUE))
list(RANK = ~ qnorm((rank(AERNASE.clin.hdac9[,var.temp], na.last = "keep") - 0.5) / sum(!is.na(AERNASE.clin.hdac9[,var.temp]))))
)
# str(UCORBIOGSAqc$Z)
cat(paste0("* renaming RANK to ", var.temp.rank,".\n"))
AERNASE.clin.hdac9[,var.temp.rank] <- NULL
names(AERNASE.clin.hdac9)[names(AERNASE.clin.hdac9) == "RANK"] <- var.temp.rank
}
# rm(var.temp, var.temp.rank)We will just visualize these transformations.
proteins_of_interest_rank_target <- c("HDAC9", proteins_of_interest_rank)
proteins_of_interest_target <- c("HDAC9", proteins_of_interest)
for(PROTEIN in proteins_of_interest_target){
cat(paste0("Plotting protein ", PROTEIN, ".\n"))
p1 <- ggpubr::gghistogram(AERNASE.clin.hdac9, PROTEIN,
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#1290D9", "#DB003F"),
add = "mean",
# rug = TRUE,
# add.params = list(color = "black", linetype = 2),
title = paste0(PROTEIN, " plaque levels"),
xlab = "",
ggtheme = theme_minimal())
print(p1)
}
for(PROTEIN in proteins_of_interest_rank_target){
cat(paste0("Plotting protein ", PROTEIN, ".\n"))
p1 <- ggpubr::gghistogram(AERNASE.clin.hdac9, PROTEIN,
# y = "..count..",
color = "white",
fill = "Gender",
palette = c("#1290D9", "#DB003F"),
add = "mean",
# rug = TRUE,
# add.params = list(color = "black", linetype = 2),
title = paste0(PROTEIN, " plaque levels"),
xlab = "inverse-normal transformation",
ggtheme = theme_minimal())
print(p1)
}
Here we calculate correlations between target(s) and 28 other cytokines. We use Spearman’s test, thus, correlations a given in rho. Please note the indications of measurement methods:
# Installation of ggcorrplot()
# --------------------------------
if(!require(devtools))
install.packages("devtools")
devtools::install_github("kassambara/ggcorrplot")
library(ggcorrplot)
# Creating matrix - inverse-rank transformation
# --------------------------------
temp <- subset(AERNASE.clin.hdac9,
select = c(proteins_of_interest_rank_target)
)
# str(AEDB.CEA.temp)
matrix.RANK <- as.matrix(temp)
rm(temp)
corr_biomarkers.rank <- round(cor(matrix.RANK,
use = "pairwise.complete.obs", #the correlation or covariance between each pair of variables is computed using all complete pairs of observations on those variables
method = "spearman"), 3)
# corr_biomarkers.rank
rename_proteins_of_interest_target <- c("HDAC9 (RNA)",
"IL2", "IL4", "IL5", "IL6", "IL8", "IL9", "IL10", "IL12",
"IL13 (L)", "IL21 (L)",
"INFG", "TNFA", "MIF (L)",
"MCP1 (L)", "MIP1a (L)", "RANTES (L)", "MIG (L)", "IP10 (L)",
"Eotaxin1 (L)", "TARC (L)", "PARC (L)", "MDC (L)",
"OPG (L)", "sICAM1 (L)", "VEGFA (E)", "TGFB (E)", "MMP2 (a)", "MMP8 (a)", "MMP9 (a)")
colnames(corr_biomarkers.rank) <- c(rename_proteins_of_interest_target)
rownames(corr_biomarkers.rank) <- c(rename_proteins_of_interest_target)
corr_biomarkers_p.rank <- ggcorrplot::cor_pmat(matrix.RANK, use = "pairwise.complete.obs", method = "spearman")
# ++++++++++++++++++++++++++++
# flattenCorrMatrix
# ++++++++++++++++++++++++++++
# cormat : matrix of the correlation coefficients
# pmat : matrix of the correlation p-values
flattenCorrMatrix <- function(cormat, pmat) {
ut <- upper.tri(cormat)
data.frame(
row = rownames(cormat)[row(cormat)[ut]],
column = rownames(cormat)[col(cormat)[ut]],
cor =(cormat)[ut],
p = pmat[ut]
)
}
corr_biomarkers.rank.df <- flattenCorrMatrix(corr_biomarkers.rank, corr_biomarkers_p.rank)
names(corr_biomarkers.rank.df)[names(corr_biomarkers.rank.df) == "row"] <- "Cytokine_X"
names(corr_biomarkers.rank.df)[names(corr_biomarkers.rank.df) == "column"] <- "CytokineY"
names(corr_biomarkers.rank.df)[names(corr_biomarkers.rank.df) == "cor"] <- "SpearmanRho"
DT::datatable(corr_biomarkers.rank.df)
fwrite(corr_biomarkers.rank.df, file = paste0(OUT_loc, "/",Today,".correlation_cytokines.txt"))# Add correlation coefficients
# --------------------------------
# argument lab = TRUE
p1 <- ggcorrplot(corr_biomarkers.rank,
method = "square",
type = "lower",
title = "Cross biomarker correlations",
show.legend = TRUE, legend.title = bquote("Spearman's"~italic(rho)),
ggtheme = ggplot2::theme_minimal, outline.color = "#FFFFFF",
show.diag = TRUE,
hc.order = FALSE,
lab = FALSE,
digits = 3,
tl.cex = 16,
# xlab = c("MCP1"),
# p.mat = corr_biomarkers_p.rank, sig.level = 0.05,
colors = c("#1290D9", "#FFFFFF", "#E55738"))
p1
ggsave(filename = paste0(PLOT_loc, "/", Today, ".correlation_cytokines.png"), plot = last_plot())
ggsave(filename = paste0(PLOT_loc, "/", Today, ".correlation_cytokines.pdf"), plot = last_plot())
rm(p1)While visually attractive we are not necessarily interested in the correlations between all the cytokines, rather of target(s)` with other cytokines only.
temp <- subset(corr_biomarkers.rank.df, Cytokine_X == "HDAC9 (RNA)" )
temp$p_log10 <- -log10(temp$p)
p_threshold <- -log10(0.05/nrow(temp))
p_threshold
p1 <- ggbarplot(temp, x = "CytokineY", y = "SpearmanRho",
fill = "CytokineY", # change fill color by cyl
# color = "white", # Set bar border colors to white
palette = uithof_color, # jco journal color palett. see ?ggpar
xlab = "Cytokine",
ylab = expression("Spearman's"~italic(rho)),
sort.val = "desc", # Sort the value in dscending order
sort.by.groups = FALSE, # Don't sort inside each group
x.text.angle = 45, # Rotate vertically x axis texts
cex = 1.25
)
ggpar(p1, legend = "bottom",
legend.title = "") +
theme(axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
axis.title.x = element_text(size = 18),
axis.title.y = element_text(size = 18))
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,"_vs_Cytokines.png"), plot = last_plot())
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,"_vs_Cytokines.pdf"), plot = last_plot())
rm(p1)Another version - probably not good.
temp <- subset(corr_biomarkers.rank.df, Cytokine_X == "HDAC9 (RNA)" )
temp$p_log10 <- -log10(temp$p)
p_threshold <- -log10(0.05/nrow(temp))
p_threshold
p1 <- ggdotchart(temp, x = "CytokineY", y = "p_log10",
color = "CytokineY", #fill = "CytokineY", # Color by groups
palette = uithof_color, # Custom color palette
xlab = "Cytokine",
ylab = expression(log[10]~"("~italic(p)~")-value"),
# ylim = c(0, 9),
sorting = "descending", # Sort value in descending order
add = "segments", # Add segments from y = 0 to dots
rotate = FALSE, # Rotate vertically
# group = "CytokineY", # Order by groups
dot.size = 16, # Large dot size
label = round(temp$SpearmanRho, digits = 3), # Add mpg values as dot labels
font.label = list(color = "white", size = 12,
vjust = 0.5)
)
ggpar(p1, legend = "",
legend.title = "") +
theme(axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
axis.title.x = element_text(size = 18),
axis.title.y = element_text(size = 18))
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,"_vs_Cytokines.dotchart.png"), plot = last_plot())
ggsave(file = paste0(PLOT_loc, "/",Today,".AERNASE.clin.hdac9.",TRAIT_OF_INTEREST,"_vs_Cytokines.dotchart.pdf"), plot = last_plot())
rm(temp, p1)In this model we correct for Age, Gender, and year of surgery.
Here we use the inverse-rank normalized data - visually this is more normally distributed.
Analysis of plaque cytokines traits as a function of plaque target(s) levels.
GLM.results <- data.frame(matrix(NA, ncol = 15, nrow = 0))
cat("Running linear regression...\n")
for (protein in 1:length(TRAITS.TARGET.RANK)) {
PROTEIN = TRAITS.TARGET.RANK[protein]
cat(paste0("\nAnalysis of ",PROTEIN,".\n"))
for (trait in 1:length(proteins_of_interest_rank)) {
TRAIT = proteins_of_interest_rank[trait]
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., PROTEIN, TRAIT, COVARIATES_M1) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
### univariate
# fit <- lm(currentDF[,PROTEIN] ~ currentDF[,TRAIT] + Age + Gender + ORdate_year, data = currentDF)
fit <- lm(currentDF[,PROTEIN] ~ currentDF[,TRAIT] + Age + Gender + ORdate_epoch, data = currentDF)
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 15, nrow = 0))
GLM.results.TEMP[1,] = GLM.CON(fit, "AEDB.CEA", PROTEIN, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
}
cat("Edit the column names...\n")
colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"T-value", "P-value", "r^2", "r^2_adj", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")
GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`T-value` <- as.numeric(GLM.results$`T-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2` <- as.numeric(GLM.results$`r^2`)
GLM.results$`r^2_adj` <- as.numeric(GLM.results$`r^2_adj`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)DT::datatable(GLM.results)
# Save the data
cat("Writing results to Excel-file...\n")
### Univariate
library(openxlsx)
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Con.Uni.",TRAIT_OF_INTEREST,"_Plaque.Cytokines_Plaques.RANK.MODEL1.xlsx"),
rowNmes = FALSE, colNames = TRUE, sheetName = "Con.Uni.PlaquePheno")
# Removing intermediates
cat("Removing intermediate files...\n")
rm(TRAIT, trait, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)In this model we correct for Age, Gender, year of surgery, Hypertension status, Diabetes status, current smoker status, lipid-lowering drugs (LLDs), antiplatelet medication, eGFR (MDRD), BMI, MedHx_CVD (combination of CAD history, stroke history, and peripheral interventions), and stenosis.
Here we use the inverse-rank normalized data - visually this is more normally distributed.
Analysis of plaque cytokines as a function of plaque target(s) levels.
GLM.results <- data.frame(matrix(NA, ncol = 15, nrow = 0))
cat("Running linear regression...\n")
for (protein in 1:length(TRAITS.TARGET.RANK)) {
PROTEIN = TRAITS.TARGET.RANK[protein]
cat(paste0("\nAnalysis of ",PROTEIN,".\n"))
for (trait in 1:length(proteins_of_interest_rank)) {
TRAIT = proteins_of_interest_rank[trait]
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., PROTEIN, TRAIT, COVARIATES_M2) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.)))
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
### univariate
# fit <- lm(currentDF[,PROTEIN] ~ currentDF[,TRAIT] + Age + Gender + ORdate_year +
# Hypertension.composite + DiabetesStatus + SmokerStatus +
# Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
# MedHx_CVD + stenose,
# data = currentDF)
fit <- lm(currentDF[,PROTEIN] ~ currentDF[,TRAIT] + Age + Gender + ORdate_epoch +
Hypertension.composite + DiabetesStatus + SmokerStatus +
Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
MedHx_CVD + stenose,
data = currentDF)
model_step <- stepAIC(fit, direction = "both", trace = FALSE)
print(model_step)
print(summary(fit))
GLM.results.TEMP <- data.frame(matrix(NA, ncol = 15, nrow = 0))
GLM.results.TEMP[1,] = GLM.CON(fit, "AEDB.CEA", PROTEIN, TRAIT, verbose = TRUE)
GLM.results = rbind(GLM.results, GLM.results.TEMP)
}
}
cat("Edit the column names...\n")
colnames(GLM.results) = c("Dataset", "Predictor", "Trait",
"Beta", "s.e.m.",
"OR", "low95CI", "up95CI",
"T-value", "P-value", "r^2", "r^2_adj", "N", "Model_N", "Perc_Miss")
cat("Correct the variable types...\n")
GLM.results$Beta <- as.numeric(GLM.results$Beta)
GLM.results$s.e.m. <- as.numeric(GLM.results$s.e.m.)
GLM.results$OR <- as.numeric(GLM.results$OR)
GLM.results$low95CI <- as.numeric(GLM.results$low95CI)
GLM.results$up95CI <- as.numeric(GLM.results$up95CI)
GLM.results$`T-value` <- as.numeric(GLM.results$`T-value`)
GLM.results$`P-value` <- as.numeric(GLM.results$`P-value`)
GLM.results$`r^2` <- as.numeric(GLM.results$`r^2`)
GLM.results$`r^2_adj` <- as.numeric(GLM.results$`r^2_adj`)
GLM.results$`N` <- as.numeric(GLM.results$`N`)
GLM.results$`Model_N` <- as.numeric(GLM.results$`Model_N`)
GLM.results$`Perc_Miss` <- as.numeric(GLM.results$`Perc_Miss`)DT::datatable(GLM.results)
# Save the data
cat("Writing results to Excel-file...\n")
### Univariate
library(openxlsx)
write.xlsx(GLM.results,
file = paste0(OUT_loc, "/",Today,".AERNASE.clin.hdac9.Con.Multi.",TRAIT_OF_INTEREST,"_Plaque.Cytokines_Plaques.RANK.MODEL2.xlsx"),
rowNames = FALSE, colNames = TRUE, sheetName = "Con.Multi.PlaquePheno")
# Removing intermediates
cat("Removing intermediate files...\n")
rm(TRAIT, trait, currentDF, GLM.results, GLM.results.TEMP, fit, model_step)Here we plot the levels of inverse-rank normal transformed target(s)
plaque levels from experiment 1 and 2 to the
Plaque vulnerability index.
library(sjlabelled)
AERNASE.clin.hdac9$yeartemp <- as.numeric(year(AERNASE.clin.hdac9$dateok))
attach(AERNASE.clin.hdac9)
AERNASE.clin.hdac9[,"ORyearGroup"] <- NA
AERNASE.clin.hdac9$ORyearGroup[yeartemp <= 2007] <- "< 2007"
AERNASE.clin.hdac9$ORyearGroup[yeartemp > 2007] <- "> 2007"
detach(AERNASE.clin.hdac9)
table(AERNASE.clin.hdac9$ORyearGroup, AERNASE.clin.hdac9$ORdate_year)# Global test
compare_means(HDAC9 ~ Plaque_Vulnerability_Index, data = AERNASE.clin.hdac9, method = "kruskal.test")
p1 <- ggpubr::ggboxplot(AERNASE.clin.hdac9,
x = "Plaque_Vulnerability_Index",
y = "HDAC9",
xlab = "Plaque vulnerability index",
ylab = "HDAC9 (normalized expression)\n",
color = "Plaque_Vulnerability_Index",
palette = "npg",
add = "jitter",
add.params = list(size = 2, jitter = 0.2)) +
stat_compare_means(label = "p.format", method = "kruskal.test") +
font("xlab", size = 17) +
font("ylab", size = 17) +
font("xy.text", size = 16) +
font("legend.title", face = "bold")
ggpar(p1, legend = "bottom", legend.title = "Plaque vulnerability index")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.PlaqueVulnerabilityIndex.pdf"), plot = last_plot())# Global test
compare_means(HDAC9 ~ Plaque_Vulnerability_Index, data = AERNASE.clin.hdac9, method = "kruskal.test")
p1 <- ggpubr::ggbarplot(AERNASE.clin.hdac9,
x = "Plaque_Vulnerability_Index",
y = "HDAC9",
xlab = "Plaque vulnerability index",
ylab = "HDAC9 (normalized expression)\n",
col = "Plaque_Vulnerability_Index",
fill = "Plaque_Vulnerability_Index",
palette = "npg",
add = "median_iqr", error.plot = "upper_errorbar") +
stat_compare_means(label = "p.format", method = "kruskal.test",
label.x = 1, label.y = 50) +
font("xlab", size = 17) +
font("ylab", size = 17) +
font("xy.text", size = 16) +
font("legend.title", face = "bold")
ggpar(p1, legend = "bottom", legend.title = "Plaque vulnerability index", ylim = c(0, 55))
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.PlaqueVulnerabilityIndex.BarPlot.median_iqr.pdf"), plot = last_plot())
compare_means(HDAC9 ~ Plaque_Vulnerability_Index, data = AERNASE.clin.hdac9, method = "kruskal.test")p1 <- ggpubr::ggbarplot(AERNASE.clin.hdac9,
x = "Plaque_Vulnerability_Index",
y = "HDAC9",
xlab = "Plaque vulnerability index",
ylab = "HDAC9 (normalized expression)\n",
col = "Plaque_Vulnerability_Index",
fill = "Plaque_Vulnerability_Index",
palette = "npg",
add = "mean_se", error.plot = "upper_errorbar") +
stat_compare_means(label = "p.format", method = "kruskal.test",
label.x = 1, label.y = 50) +
font("xlab", size = 17) +
font("ylab", size = 17) +
font("xy.text", size = 16) +
font("legend.title", face = "bold")
ggpar(p1, legend = "bottom", legend.title = "Plaque vulnerability index", ylim = c(0, 55))
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.PlaqueVulnerabilityIndex.BarPlot.means_se.pdf"), plot = last_plot())compare_means(HDAC9 ~ Plaque_Vulnerability_Index, data = subset(AERNASE.clin.hdac9, HDAC9 <100), method = "kruskal.test")
p1 <- ggpubr::ggboxplot(subset(AERNASE.clin.hdac9, HDAC9 <100) ,
x = "Plaque_Vulnerability_Index",
y = "HDAC9",
xlab = "Plaque vulnerability index",
ylab = "HDAC9 (normalized expression)\noutliers above 100 are removed",
col = "Plaque_Vulnerability_Index",
fill = "Plaque_Vulnerability_Index",
palette = "npg",
add = "boxplot", error.plot = "crossbar") +
stat_compare_means(label = "p.format", method = "kruskal.test",
label.x = 1, label.y = 50) +
font("xlab", size = 17) +
font("ylab", size = 17) +
font("xy.text", size = 16) +
font("legend.title", face = "bold")
ggpar(p1, legend = "bottom", legend.title = "Plaque vulnerability index")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.PlaqueVulnerabilityIndex.Boxplot.outlier_above_100_removed.pdf"), plot = last_plot())compare_means(HDAC9 ~ Plaque_Vulnerability_Index, data = AERNASE.clin.hdac9, method = "kruskal.test")
p1 <- ggpubr::ggboxplot(AERNASE.clin.hdac9,
x = "Plaque_Vulnerability_Index",
y = "HDAC9",
facet.by = "Plaque_Vulnerability_Index",
xlab = "Plaque vulnerability index",
ylab = "HDAC9 (normalized expression)\n",
color = "Plaque_Vulnerability_Index",
palette = "npg",
add = "jitter",
add.params = list(size = 2, jitter = 0.2)) +
stat_compare_means(label = "p.format", method = "kruskal.test") +
font("xlab", size = 17) +
font("ylab", size = 17) +
font("xy.text", size = 16) +
font("legend.title", face = "bold")
ggpar(p1, legend = "bottom", legend.title = "Plaque vulnerability index")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.FacetByPlaqueVulnerabilityIndex.pdf"), plot = last_plot())compare_means(HDAC9 ~ Plaque_Vulnerability_Index, group.by = "Gender", data = AERNASE.clin.hdac9, method = "kruskal.test")
p2 <- ggpubr::ggboxplot(AERNASE.clin.hdac9,
x = "Plaque_Vulnerability_Index",
y = "HDAC9",
xlab = "Plaque vulnerability index by gender",
ylab = "HDAC9 (normalized expression)\n",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")
ggpar(p2, legend = "bottom", legend.title = "Plaque vulnerability index")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.PlaqueVulnerabilityIndex.byGender.pdf"), plot = last_plot())
compare_means(HDAC9 ~ Plaque_Vulnerability_Index, data = AERNASE.clin.hdac9, method = "kruskal.test")
p5 <- ggpubr::ggboxplot(AERNASE.clin.hdac9,
x = "Plaque_Vulnerability_Index",
y = "HDAC9",
xlab = "Plaque vulnerability index",
ylab = "HDAC9 (normalized expression)\n",
color = "Plaque_Vulnerability_Index",
palette = "npg",
facet.by = "ORyearGroup",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")
ggpar(p5, legend = "bottom", legend.title = "Plaque vulnerability index")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.PlaqueVulnerabilityIndex_Facet_byYear.pdf"), plot = last_plot())compare_means(HDAC9 ~ Plaque_Vulnerability_Index, group.by = "Gender", data = AERNASE.clin.hdac9, method = "kruskal.test")
p6 <- ggpubr::ggboxplot(AERNASE.clin.hdac9,
x = "Plaque_Vulnerability_Index",
y = "HDAC9",
xlab = "Plaque vulnerability index",
ylab = "HDAC9 (normalized expression)\n",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
facet.by = "ORyearGroup",
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")
ggpar(p6, legend = "bottom", legend.title = "Plaque vulnerability index")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.PlaqueVulnerabilityIndex_Facet_byYear.byGender.pdf"), plot = last_plot())In this model we correct for Age, Gender, and year of surgery.
Here we use the inverse-rank normalized data - visually this is more normally distributed.
Analysis of the plaque vulnerability indez as a function of plaque target(s) levels.
TRAITS.TARGET.RANK.extra = c("HDAC9")
GLM.results <- data.frame(matrix(NA, ncol = 16, nrow = 0))
for (protein in 1:length(TRAITS.TARGET.RANK.extra)) {
PROTEIN = TRAITS.TARGET.RANK.extra[protein]
cat(paste0("\nAnalysis of ",PROTEIN,".\n"))
TRAIT = "Plaque_Vulnerability_Index"
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., PROTEIN, TRAIT, COVARIATES_M1) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.))) %>%
droplevels(.)
# fix numeric OR year
# currentDF$ORdate_year <- as.numeric(currentDF$ORdate_year)
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
# print(class(currentDF[,TRAIT]))
# table(currentDF$ORdate_year)
### univariate
# + Hypertension.composite + DiabetesStatus + SmokerCurrent +
# Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI +
# CAD_history + Stroke_history + Peripheral.interv + stenose
# fit <- polr(currentDF[,TRAIT] ~ currentDF[,PROTEIN] + Age + Gender + ORdate_year,
# data = currentDF,
# Hess = TRUE)
fit <- polr(currentDF[,TRAIT] ~ currentDF[,PROTEIN] + Age + Gender + ORdate_epoch,
data = currentDF,
Hess = TRUE)
print(summary(fit))
## store table
(ctable <- coef(summary(fit)))
## calculate and store p values
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
## combined table
print((ctable <- cbind(ctable, "p value" = p)))
}In this model we correct for Age, Gender, Hypertension status, Diabetes status, current smoker status, lipid-lowering drugs (LLDs), antiplatelet medication, eGFR (MDRD), BMI, MedHx_CVD (combination of CAD history, stroke history, and peripheral interventions), and stenosis..
for (protein in 1:length(TRAITS.TARGET.RANK.extra)) {
PROTEIN = TRAITS.TARGET.RANK.extra[protein]
cat(paste0("\nAnalysis of ",PROTEIN,".\n"))
TRAIT = "Plaque_Vulnerability_Index"
cat(paste0("\n- processing ",TRAIT,"\n\n"))
currentDF <- as.data.frame(AERNASE.clin.hdac9 %>%
dplyr::select(., PROTEIN, TRAIT, COVARIATES_M2) %>%
filter(complete.cases(.))) %>%
filter_if(~is.numeric(.), all_vars(!is.infinite(.))) %>%
droplevels(.)
# fix numeric OR year
# currentDF$ORdate_year <- as.numeric(currentDF$ORdate_year)
# for debug
# print(DT::datatable(currentDF))
# print(nrow(currentDF))
# print(str(currentDF))
# print(class(currentDF[,TRAIT]))
### univariate
# fit <- polr(as.factor(currentDF[,TRAIT]) ~ currentDF[,PROTEIN] + Age + Gender + ORdate_year + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose,
# data = currentDF,
# Hess = TRUE)
fit <- polr(as.factor(currentDF[,TRAIT]) ~ currentDF[,PROTEIN] + Age + Gender + ORdate_epoch + Hypertension.composite + DiabetesStatus + SmokerStatus + Med.Statin.LLD + Med.all.antiplatelet + GFR_MDRD + BMI + MedHx_CVD + stenose,
data = currentDF,
Hess = TRUE)
print(summary(fit))
## store table
(ctable <- coef(summary(fit)))
## calculate and store p values
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
## combined table
print((ctable <- cbind(ctable, "p value" = p)))
}# Global test
compare_means(HDAC9 ~ Fat.bin_10, data = AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_10)), method = "kruskal.test")
p1 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_10)),
x = "Fat.bin_10",
y = "HDAC9",
xlab = "Fat <10% vs >10%",
ylab = "HDAC9 (normalized expression)\n",
color = "Fat.bin_10",
palette = "npg",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")
ggpar(p1, legend = "bottom", legend.title = "Fat <10% vs >10%")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.Fat.bin_10.pdf"), plot = last_plot())compare_means(HDAC9 ~ Fat.bin_10, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_10)), method = "kruskal.test")
p2 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_10)),
x = "Fat.bin_10",
y = "HDAC9",
xlab = "Fat <10% vs >10% by gender",
ylab = "HDAC9 (normalized expression)\n",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")
ggpar(p2, legend = "bottom", legend.title = "Fat <10% vs >10% by gender")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.Fat.bin_10.byGender.pdf"), plot = last_plot())compare_means(HDAC9 ~ Fat.bin_10, data = AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_10)), method = "kruskal.test")
p5 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_10)),
x = "Fat.bin_10",
y = "HDAC9",
xlab = "Fat <10% vs >10% by year of surgery",
ylab = "HDAC9 (normalized expression)\n",
color = "Fat.bin_10",
palette = "npg",
facet.by = "ORyearGroup",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")
ggpar(p5, legend = "bottom", legend.title = "Fat <10% vs >10% by year of surgery")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.Fat.bin_10_Facet_byYear.pdf"), plot = last_plot())compare_means(HDAC9 ~ Fat.bin_10, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_10)), method = "kruskal.test")
p6 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_10)),
x = "Fat.bin_10",
y = "HDAC9",
xlab = "Fat <10% vs >10% by year of surgery and gender",
ylab = "HDAC9 (normalized expression)\n",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
facet.by = "ORyearGroup",
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")
ggpar(p6, legend = "bottom", legend.title = "Fat <10% vs >10% by year of surgery and gender")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.Fat.bin_10_Facet_byYear.byGender.pdf"), plot = last_plot())# Global test
compare_means(HDAC9 ~ Fat.bin_40, data = AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_40)), method = "kruskal.test")
p1 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_40)),
x = "Fat.bin_40",
y = "HDAC9",
xlab = "Fat <40% vs >40%",
ylab = "HDAC9 (normalized expression)\n",
color = "Fat.bin_40",
palette = "npg",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")
ggpar(p1, legend = "bottom", legend.title = "Fat <40% vs >40%")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.Fat.bin_40.pdf"), plot = last_plot())compare_means(HDAC9 ~ Fat.bin_40, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_40)), method = "kruskal.test")
p2 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_40)),
x = "Fat.bin_40",
y = "HDAC9",
xlab = "Fat <40% vs >40% by gender",
ylab = "HDAC9 (normalized expression)\n",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")
ggpar(p2, legend = "bottom", legend.title = "Fat <40% vs >40% by gender")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.Fat.bin_40.byGender.pdf"), plot = last_plot())compare_means(HDAC9 ~ Fat.bin_40, data = AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_40)), method = "kruskal.test")
p5 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_40)),
x = "Fat.bin_40",
y = "HDAC9",
xlab = "Fat <40% vs >40% by year of surgery",
ylab = "HDAC9 (normalized expression)\n",
color = "Fat.bin_40",
palette = "npg",
facet.by = "ORyearGroup",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")
ggpar(p5, legend = "bottom", legend.title = "Fat <40% vs >40% by year of surgery")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.Fat.bin_40_Facet_byYear.pdf"), plot = last_plot())compare_means(HDAC9 ~ Fat.bin_40, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_40)), method = "kruskal.test")
p6 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Fat.bin_40)),
x = "Fat.bin_40",
y = "HDAC9",
xlab = "Fat <40% vs >40% by year of surgery and gender",
ylab = "HDAC9 (normalized expression)\n",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
facet.by = "ORyearGroup",
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")
ggpar(p6, legend = "bottom", legend.title = "Fat <40% vs >40% by year of surgery and gender")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.Fat.bin_40_Facet_byYear.byGender.pdf"), plot = last_plot())# Global test
compare_means(HDAC9 ~ IPH.bin, data = AERNASE.clin.hdac9 %>% filter(!is.na(IPH.bin)), method = "kruskal.test")
p1 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(IPH.bin)),
x = "IPH.bin",
y = "HDAC9",
xlab = "Intraplaque hemorrhage (no vs. yes)",
ylab = "HDAC9 (normalized expression)\n",
color = "IPH.bin",
palette = "npg",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")
ggpar(p1, legend = "bottom", legend.title = "IPH")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.IPH.bin.pdf"), plot = last_plot())compare_means(HDAC9 ~ IPH.bin, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(IPH.bin)), method = "kruskal.test")
p2 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(IPH.bin)),
x = "IPH.bin",
y = "HDAC9",
xlab = "Intraplaque hemorrhage (no vs. yes) by gender",
ylab = "HDAC9 (normalized expression)\n",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")
ggpar(p2, legend = "bottom", legend.title = "IPH by gender")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.IPH.bin.byGender.pdf"), plot = last_plot())compare_means(HDAC9 ~ IPH.bin, data = AERNASE.clin.hdac9 %>% filter(!is.na(IPH.bin)), method = "kruskal.test")
p5 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(IPH.bin)),
x = "IPH.bin",
y = "HDAC9",
xlab = "Intraplaque hemorrhage (no vs. yes) by year of surgery",
ylab = "HDAC9 (normalized expression)\n",
color = "IPH.bin",
palette = "npg",
facet.by = "ORyearGroup",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")
ggpar(p5, legend = "bottom", legend.title = "IPH by year of surgery")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.IPH.bin_Facet_byYear.pdf"), plot = last_plot())compare_means(HDAC9 ~ IPH.bin, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(IPH.bin)), method = "kruskal.test")
p6 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(IPH.bin)),
x = "IPH.bin",
y = "HDAC9",
xlab = "Intraplaque hemorrhage (no vs. yes) by year of surgery and gender",
ylab = "HDAC9 (normalized expression)\n",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
facet.by = "ORyearGroup",
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")
ggpar(p6, legend = "bottom", legend.title = "IPH by year of surgery and gender")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.IPH.bin_Facet_byYear.byGender.pdf"), plot = last_plot())# Global test
compare_means(HDAC9 ~ Calc.bin, data = AERNASE.clin.hdac9 %>% filter(!is.na(Calc.bin)), method = "kruskal.test")
p1 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Calc.bin)),
x = "Calc.bin",
y = "HDAC9",
xlab = "Calcification (no/minor vs. moderate/heavy)",
ylab = "HDAC9 (normalized expression)\n",
color = "Calc.bin",
palette = "npg",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")
ggpar(p1, legend = "bottom", legend.title = "Calcification")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.Calc.bin.pdf"), plot = last_plot())compare_means(HDAC9 ~ Calc.bin, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(Calc.bin)), method = "kruskal.test")
p2 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Calc.bin)),
x = "Calc.bin",
y = "HDAC9",
xlab = "Calcification (no/minor vs. moderate/heavy) by gender",
ylab = "HDAC9 (normalized expression)\n",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")
ggpar(p2, legend = "bottom", legend.title = "Calcification by gender")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.Calc.bin.byGender.pdf"), plot = last_plot())compare_means(HDAC9 ~ Calc.bin, data = AERNASE.clin.hdac9 %>% filter(!is.na(Calc.bin)), method = "kruskal.test")
p5 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Calc.bin)),
x = "Calc.bin",
y = "HDAC9",
xlab = "Calcification (no/minor vs. moderate/heavy) by year of surgery",
ylab = "HDAC9 (normalized expression)\n",
color = "Calc.bin",
palette = "npg",
facet.by = "ORyearGroup",
add = "jitter") +
stat_compare_means(label = "p.format", method = "kruskal.test")
ggpar(p5, legend = "bottom", legend.title = "Calcification by year of surgery")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.Calc.bin_Facet_byYear.pdf"), plot = last_plot())compare_means(HDAC9 ~ Calc.bin, group.by = "Gender", data = AERNASE.clin.hdac9 %>% filter(!is.na(Calc.bin)), method = "kruskal.test")
p6 <- ggpubr::ggboxplot(AERNASE.clin.hdac9 %>% filter(!is.na(Calc.bin)),
x = "Calc.bin",
y = "HDAC9",
xlab = "Calcification (no/minor vs. moderate/heavy) by year of surgery and gender",
ylab = "HDAC9 (normalized expression)\n",
color = "Gender",
palette = c("#D5267B", "#1290D9"),
facet.by = "ORyearGroup",
add = "jitter") +
stat_compare_means(aes(group = Gender), label = "p.format", method = "kruskal.test")
ggpar(p6, legend = "bottom", legend.title = "Calcification by year of surgery and gender")
ggsave(filename = paste0(PLOT_loc, "/", Today, ".",TRAIT_OF_INTEREST,".plaque.Calc.bin_Facet_byYear.byGender.pdf"), plot = last_plot())Version: v1.0.5
Last update: 2023-05-31
Written by: Sander W. van der Laan (s.w.vanderlaan-2[at]umcutrecht.nl).
Description: Script to analyse HDAC9 from the Ather-Express Biobank Study.
Minimum requirements: R version 3.5.2 (2018-12-20) -- 'Eggshell Igloo', macOS Mojave (10.14.2).
**MoSCoW To-Do List**
The things we Must, Should, Could, and Would have given the time we have.
_M_
_S_
_C_
_W_
**Changes log**
* v1.0.5 Fixed forest plot and alternative boxplot for symptoms.
* v1.0.4 Made histogram of PVI. Exported HDAC9 and PVI data.
* v1.0.3 Small adaptations to PVI-plots.
* v1.0.2 Changed the PVI-plot.
* v1.0.1 Added figures on fat in the plaque.
* v1.0.0 Inital version.
sessionInfo()save.image(paste0(PROJECT_loc, "/",Today,".",PROJECTNAME,".bulkRNAseq.additional_figures.RData"))| © 1979-2023 Sander W. van der Laan | s.w.vanderlaan[at]gmail.com | vanderlaan.science. | |